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ABSTRACT 

We report on a problem found in Mercury, a hybrid symplectic integrator used for 
dynamical problems in Astronomy. The variable that keeps track of bodies' statuses 
is uninitialised, which can result in bodies disappearing from simulations in a non- 
physical manner. Some FORTRAN compilers implicitly initialise variables, preventing 
simulations from having this problem. With others compilers, simulations with a suit- 
ably large maximum number of bodies parameter value are also unaffected. Otherwise, 
the problem manifests at the first event after the integrator is started, whether from 
scratch or continuing a previously stopped simulation. Although the problem does not 
manifest in some conditions, explicitly initialising the variable solves the problem in 
a permanent and unconditional manner. 
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1 INTRODUCTION 



Mercury (|Chambersl Il999l ) is a general-purpose software 
package, written in FORTRAN77, for doing N-body integra- 
tions to investigate dynamical problems in Astronomy. In a 
set of simulations to study accretion dy namics using Mer- 
cury |de Souza Torres fc Winter! 120081 ). it was found that 
the results contained fewer planets than expected. Analysis 
of output files and creation of scripts and movies to carefully 
follow the processes showed discontinuities in the number of 
embryos during the integrations. No events for these bodies 
were being registered in the output files and their disappear- 
ance was non-physical. 



Table 1. Valid values for the variable STAT in Mercury integra- 
tor. 



STAT value 


Body status 





Alive 


-2 


Collided 


-3 


Ejected 



values in its respective STAT elements, causing non-physical 
disappearance of some bodies during the integrations. 



2 THE PROBLEM 

In mercury, an integer variable array called STAT keeps 
track of the statuses of particles: whether they are alive, 
have been involved in a collision or have been ejected from 
the system. At the beginning of an old or new integration, 
all the bodies must have its STAT equals to zero (Table [TJ. 
When a live body is collided or ejected, its STAT value is 
changed accordingly and the body is removed by the sub- 
routine called mxx_elim. 

In the code, the STAT array is not initialised, meaning 
its elements may access random values. If these values are 
negative, not only bodies that have been involved in an event 
are going to be removed, but also those with invalid negative 
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2.1 Characterising the problem 

Not all simulations are affected by this problem. Some FOR- 
TRAN compilers (e.g. ifort) i mplicitly initialise variables o n 
their declaration statement (jChivers fc Sleightholmel l2006h 
Using one of them, an integer array would have all its ele- 
ments set to zero in its declaration, even without an explicit 
initialisation in the source code. Others compilers (e.g. g77 
and gfortran) begin initialising the early elements of an inte- 
ger array with zero when the array is beyond a certain size; 
this size will be compiler and machine environment depen- 
dent. With these compilers, the results will be affected by 
the non-initialisation of the STAT array if uninitialised ele- 
ments are used; it is this we are about to investigate further. 

The STAT array has its size initially defined by the pa- 
rameter nmax, a maximum number of bodies set by the user 
in the configuration file mercury. INC. During execution of 
mercury, the number of STAT elements used is equivalent 
to the actual number of bodies (nbod) in the simulation. If 
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STAT had been initialised then nmax would need merely to 
be equal to or greater than nbod. However, when using a 
compiler such as g77 or gfortran, problems will be encoun- 
tered unless nmax is greater than nbod by several hundreds. 

Using a g77 compiler, simulations with 25 different ini- 
tial conditions (Tabic [2 1), including the example given by 
Mercury's author (Chambers 2001J), were ran one or more 
times with different values of nmax. An audit was then con- 
ducted of the number of bodies and a pattern emerged: sim- 
ulations with a value of nmax considerably higher than the 
number of bodies did not have problems; bodies disappeared 
from simulations with a comparatively small nmax value. 

A set of 12 simulations, with different numbers of bod- 
ies, were used to find the minimum limit of nmax required 
for a problem-free execution (depicted with a * in Table [5] 
Two different environments were used, a Cygwin with GNU 
Bash, version 3.2.39(19)-release, and a Debian 4.0 with ker- 
nel 2.6.18-5-amd64. For each execution, a short integration 
time (about 100 years) was used. The value of nmax was 
varied in each simulation until the lower limit required to 
avoid losing bodies in a non-physical manner was found. 
The lower limit was found to be proportional to the num- 
ber of bodies (Figure [TJ. As the number of bodies tends to 
zero, the lower limit of nmax tends toward a non-zero value 
around 500-700, the exact value depending on the comput- 
ing environment. 

The best fit for the points in the Cygwin environment 

is: 

f(x) = 1.00341a; + 542.708 (1) 
and in the Debian environment: 

g(x) = 1.00009a; + 637.329 (2) 

where x is the total number of bodies. These functions are 
basically the number of bodies plus an offset value. These 
tests offer a basic way to verify if old or recent simulations 
could be affected for the problem in STAT variable's values. 
We advise that any simulation with the value of the param- 
eter nmax close to these limits to have its results checked, 
if the compiler used was g77, gfortran or one with similar 
features. 

A simple test program was run, in which a variable array 
was declared but not initialised and its values output. The 
array values were mostly large positive, large negative and 
zero. Varying the size of the array, it was found that beyond 
a limiting size the values at the beginning of the array were 0. 
Figure[2]shows this relation when using the g77 compiler in a 
Debian environment; when the number of elements exceeded 
696, zeros began appearing at the beginning of the array. 
Therefore, if nmax is greater than the sum of this minimum 
limit plus the number of bodies then no bodies will disappear 
non-physically. 

The test was repeated with the gfortran and ifort com- 
pilers. With gfortran the minimum limit was found to be 
628. With ifort the minimum value was found to be as the 
compiler initialises all array elements to 0. 



3 SOLUTION 

Excepting when the used FORTRAN compiler is one that 
makes unconditional implicitly initialisation, a variable not 



explicitly initialised can have unpredictable behavior. In or- 
der to solve this problem, the variable array STAT must 
be initialised with zero in some point inside the file Mer- 
cury6-2.for, before the commands' block for the main calcu- 
lations. We propose this initialisation to be done in the sub- 
routine MIOJN, before (line 6046) or after (line 6111) the 
block of commands "Check for attempts to do incompatible 
things". Adding the command outside of the if statement 
(lines 5911-6045) for new and old integrations will guarantee 
it is going to be executed for both of them. The initialisation 
could be done with the three lines: 

do j=2, nbod 

STAT(j) = 
end do 

With this command being executed each time an inte- 
gration starts afresh or from dump files, the variable STAT 
will not receive random values, independently of the value 
of the parameter nmax (it must still respect the basic rule: 
nmax ^ nbod), the machine environment or the FORTRAN 
compiler. 

The initialisation could be done also in other points of 
the code, but one must be sure it is being done before the 
main calculations for both old and new integrations. 



4 CONCLUSIONS 

We repeated all the tests showed in Section 12.11 with a 
corrected version of Mercury, initialising the variable ar- 
ray STAT at the end of the subroutine MIOJN in the 
source code. No discontinuities were seen in any results, 
independently of conditions. The new version is now be- 
ing used in 40 new simulations of accretion dynamics. We 
believe this small change can improve the program mak- 
ing it more reliable for any type of N-body problem sim- 
ulations. The corrected version can be downloaded from 
http : //www. astro .keele . ac.uk/~dra/mercury/ 
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Table 2. Tests made to find discontinuities problem in Mercury's results. a 



Big Bodies 


Small Bodies 


Giant Planets 


Collisions 


Central 


Style 


Algorithm 


Times 





1 





no 


Sun 


Cartesian 


Hybrid 


x 1 


*1 








no 


Sun 


Cartesian 


Hybrid 


x 1 


*1 


1 





no 


Sun 


Cartesian 


MVS 


X 1 


9 


204 


1 


no 


Sun 


Cartesian 


BS 


x 3 


9 


204 


4 


yes 


Sun 


Cartesian 


Hybrid 


x 1 


9 





4 


no 


Sun 


Cartesian 


Hybrid 


x 1 


*11 


1 





no 


Jupiter 


Asteroidal 


BS 


x 2 


11 


1 





yes 


Jupiter 


Asteroidal 


Hybrid 


x 2 


*11 


4 





yes 


Jupiter 


Asteroidal 


Hybrid 


x 1 


*11 


6 





no 


Jupiter 


Asteroidal 


Hybrid 


x 1 


11 


6 


<1 


yes 


Sun 


Asteroidal 


BS2 


x 2 


*14 


2 


<1 


yes 


Sun 


Asteroidal 


BS2 


X 1 


18 


200 


4 


yes 


Sun 


Asteroidal 


MVS 


X 1 


68 


204 


1 


yes 


Sun 


Asteroidal 


Hybrid 


x 2 


69 


204 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 3 


72 


204 


1 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 


73 


204 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 3 


88 


204 


1 


yes 


Sun 


Asteroidal 


Hybrid 


x 3 


*88 


1525 


1 


yes 


Sun 


Asteroidal 


Hybrid 


x 2 


89 


204 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 9 


*89 


500 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 


*89 


800 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 


*89 


1000 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 


*89 


1300 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 


*89 


1450 


2 


yes 


Sun 


Asteroidal 


Hybrid 


x 1 



a Columns are: Number of big bodies, number of small bodies, number of giant planets, if accepts collisions or not, central 
body, style of data input, integrator algorithm and times the same instance of the simulation were simulated with small 
difference in bodies' positions. Simulations for Figure[T]are shown by *. 
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Figure 1. Minimum value of NMAX required for a problem-free execution as a function of number of bodies. The value depends somewhat 
on the environment; two were tested: Cygwin with GNU Bash, version 3.2.39(19)-release, and Debian 4.0 with kernel 2.6.18-5-amd64. 
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Figure 2. The number of elements that are and are not initialised to zero as a function of the number of elements, when using the g77 
compiler under a Debian environment. As the array size increases (abscissa), the number of uninitialised elements increases steadily, with 
a smaller number of initialised elements spread throughout the array. When the number of elements passes 696, the number uninitialised 
remains constant and the number initialised increases linearly. Beyond this value, initialised elements are added at the beginning of the 
array, pushing the uninitialised elements to the end. The result is similar when using gfortran but the turnoff occurs at 628 elements. If 
Mercury has been compiled with such a compiler then the NMAX value would have to have been greater than the turnoff value plus the 
number of bodies in order for bodies not to have disappeared in a non-physical manner. 



